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Abstract.  We  consider  constraint  preserving  multidimensional  evolution  equations. 
A  prototypical  example  is  provided  by  the  magnetic  induction  equation  of  plasma 
physics.  The  constraint  of  interest  is  the  divergence  of  the  magnetic  field.  We  de¬ 
sign  finite  volume  schemes  which  approximate  these  equations  in  a  stable  manner  and 
preserve  a  discrete  version  of  the  constraint.  The  schemes  are  based  on  reformulat¬ 
ing  standard  edge  centered  finite  volume  fluxes  in  terms  of  vertex  centered  potentials. 
The  potential-based  approach  provides  a  general  framework  for  faithful  discretizations 
of  constraint  transport  and  we  apply  it  to  both  divergence  preserving  as  well  as  curl 
preserving  equations.  We  present  benchmark  numerical  tests  which  confirm  that  our 
potential-based  schemes  achieve  high  resolution,  while  being  constraint  preserving. 

AMS  subject  classifications:  65M06,35L65 

Key  words:  multidimensional  evolution  equations,  magnetohydrodynamics,  finite-difference,  fi¬ 
nite  volume  schemes,  constraint  transport,  potential-based  schemes. 


1  Introduction 

We  are  concerned  with  evolution  equations  of  the  form 

u*  +  L(ax,f(x,f,u))  =  0,  v(x,f)e  rxE+/  (1.1) 

where  u(x,f) :  x  R+  i— >  Rm  is  the  unknown,  f :  X  is  a  nonlinear  flux  function  and 

L :  X  i— >  Y  is  a  differential  operator  acting  on  the  Sobolev  space  X.  We  assume  there  exists 
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another  differential  operator  M:Y\-^Z  such  that  ML(f(v,v))=0  for  all  v  E  X.  Applying 
the  operator  M  to  both  sides  of  (1.1),  we  obtain 


(Mu)*  =  0.  (1.2) 

Hence,  solutions  of  (1.1)  satisfy  an  additional  constraint  which  enforces  them  to  lie  on  a 
sub-manifold  of  the  space  X. 

The  above  framework  is  generic  to  a  large  class  of  evolution  equations  involving  in¬ 
trinsic  constraints.  We  mention  three  prototype  examples.  As  a  first  example,  consider  the 
curl  advection: 

Uf  +  curl(f(x,f,u))  =0,  (x,t)  eR nxR+.  (1.3) 

This  equation  is  an  example  for  (1.1),(1.2)  with  the  differential  operators  L  =  curl  and 
M  =  div.  Hence,  solutions  of  (1.3)  satisfy  the  additional  divergence  constraint: 

div(u)*  =  0.  (1.4) 

A  specific  example  for  (1.3)  is  the  magnetic  induction  equation  of  plasma  physics.  Under 

the  assumptions  of  zero  resistivity,  the  magnetic  field  u,  evolving  under  the  influence  of 
a  given  velocity  v,  satisfies  the  following  form  of  the  Maxwell's  equations  [23]: 

u^  +  curl(u  x  v)  =0,  (x,t)e RnxR+.  (1.5) 

The  fact  that  magnetic  monopoles  have  not  been  observed  in  nature  implies  that 

div(u(x,0))  =  0.  (1.6) 

As  a  consequence  of  the  divergence  constraint  (1.4),  the  solutions  of  (1.5)  remain  diver¬ 
gence  free.  The  magnetic  induction  equation  (1.5)  is  a  sub-model  for  the  equations  of 
ideal  Magnetohydrodynamics  (MHD)  [11]. 

Adding  magnetic  resistivity  to  the  model  leads  to  the  viscous  magnetic  induction 
equations: 

+  curl(u  x  v)  =  —  <7(curl(curlu)),  (x,t)  E  Rn  x  ]R+.  (1.7) 

The  parameter  a  is  the  resistivity  co-efficient  of  the  medium.  Solutions  of  (1.7)  also  satisfy 
the  divergence  constraint  (1.4). 

A  second  example  for  (1.1),(1.2)  is  the  grad  advection: 

Wf  +  grad(f(x,f,w))  =0,  (x,t)  GRnxR+.  (1.8) 

The  differential  operators  of  interest  are  L  —  grad  and  M  =  curl  and  solutions  of  (1.8) 
satisfy  the  additional  constraint. 


curl(w)^  =  0. 
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A  related  example  is  the  system  of  two  dimensional  linear  wave  equations  [21]: 

Pt  +  CUX-\~CVy  —  0, 

<  Ut  +  cpx  =  0,  (1.9) 

KVt  +  CPy  —  0. 

The  corresponding  differential  operators  L  and  M  are 

/0  c  0\  /0  0  c\ 

L=  c  0  0  \dx  +  0  0  0  \  dy, 

\0  0  0/  \c  0  0/ 

and  M=  (0,0,1  )3X+  (0,-1, 0)dy  picks  the  vorticity  of  the  flow,  co\—vx—uy. 

A  third  -  considerably  more  complicated  example  of  evolution  equations  with  non¬ 
linear  constraints  is  provided  by  the  Einstein  equations  [10]. 

Standard  finite  volume/ finite  difference  numerical  schemes  for  approximating  (1.1),(1.2) 
may  not  necessarily  treat  the  constraint  properly  and  may  fail  to  be  stable  [11].  Hence, 
suitable  constraint  preserving  schemes  need  to  be  devised  for  robust  approximation  of 

(1.1) ,(1.2).  Design  of  efficient  numerical  schemes  for  the  constrained  evolution  equations 

(1.1) ,(1.2)  is  a  highly  active  research  area,  consult  [1,9,21,24,31]  and  references  therein. 
We  mention  below  three  main  methods  available  in  the  literature  for  handling  constraint 
transport  equations,  where  most  of  the  attention  is  devoted  to  applications  to  the  ideal 
MHD  and  the  magnetic  induction  equations  (1.5). 

1.1  Pro  j  ection  method. 

This  method  [5-7]  is  based  on  the  Hodge  decomposition  of  the  solution  u  of  (1.5).  The 
update  un  at  each  time  step  may  not  be  divergence  free  and  is  corrected  by  the  decompo¬ 
sition,  un  —  VY + curlO.  Applying  the  divergence  operator  to  the  Hodge  decomposition 
leads  to  the  elliptic  equation: 

-AY  =  div(un). 

The  corrected  field  u*  =  un  —  V  Y  is  divergence  free.  This  method  can  be  very  expensive 
computationally  as  an  elliptic  equation  has  to  be  solved  at  every  time  step,  augmented 
with  proper  set  of  boundary  conditions,  e.g.  [31]. 

1.2  Source  terms. 

Adding  a  source  term  proportional  to  the  divergence  in  (1.5)  results  in 

+  curl(u  x  v)  =  —  vdiv(u).  (1.10) 

Applying  the  divergence  to  both  sides,  we  obtain 

(div(u))f+div(vdivu)  =0. 
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Hence,  any  potential  divergence  errors  are  transported  away  from  the  computational  do¬ 
main  by  the  flow.  Furthermore,  the  form  (1.10)  is  symmetrizable  [13].  This  procedure  for 
"cleaning"  the  divergence  was  introduced  in  [22,23].  Recent  papers  [11,12]  have  demon¬ 
strated  that  the  source  term  in  (1.10)  needs  to  be  discretized  in  a  very  careful  manner 
for  numerical  stability.  Another  problem  with  this  approach  lies  in  the  non-conservative 
form  of  (1.10).  Hence,  numerical  schemes  based  on  this  approach  may  result  in  wrong 
shock  speeds  [31]. 

1.3  Design  of  special  divergence  operators/staggering. 

This  popular  method  consists  of  staggering  the  discretizations  of  the  velocity  and  mag¬ 
netic  fields  in  (1.5).  A  wide  variety  of  strategies  for  staggering  the  meshes  has  been  pro¬ 
posed  [2-4,8,9,24,31]  and  references  therein.  The  presence  of  different  sets  of  meshes 
leads  to  problems  when  the  staggered  schemes  are  parallelized.  Unstaggered  variants  of 
this  approach  have  also  been  proposed  in  [1,29,30].  The  approach  suggested  in  [30]  is  of 
particular  relevance  for  this  paper;  the  authors  suggest  an  unstaggered  method,  based  on 
upwindedflux  distributions  in  each  cell,  resulting  in  a  scheme  which  preserves  a  particular 
discrete  form  of  divergence. 

The  above  examples  leave  room  for  designing  other  constraint  preserving  schemes  that 
are  easy  to  implement  and  computationally  robust.  In  this  paper  we  propose  a  new  ap¬ 
proach  for  designing  such  schemes.  Our  starting  point  is  the  genuinely  multi-dimensional 
structure  of  equation  (1.1)  complemented  with  the  constraint  (1.2).  This  is  in  contrast  to 
standard  finite  volume  schemes  based  on  locally  one  dimensional  edge  centered  fluxes 
[18, 28],  which  do  not  incorporate  any  explicit  information  in  the  transverse  direction. 
Here,  we  introduce  a  new  approach  to  modify  standard  finite  volume  schemes  which 
incorporates  genuinely  multi-dimensional  information,  resulting  in  a  new  family  of  con¬ 
straint  preserving  schemes. 

To  this  end,  we  propose  the  construction  of  vertex-centered  numerical  potentials  which 
serve  as  the  building  blocks  of  our  constraint  preserving  schemes.  Written  in  terms  of 
these  numerical  potentials,  the  proposed  schemes  are  genuinely  multi-dimensional  and 
preserve  a  discrete  form  of  the  constraint  (1.2).  The  framework  is  very  general,  easy  to 
code  and  allow  for  any  consistent  numerical  flux  to  be  used  as  a  building  block  for  the 
construction  of  numerical  potentials.  No  additional  upwinding  is  necessary  for  numeri¬ 
cal  stability.  Our  new,  so-called  potential-based  approach,  is  demonstrated  in  the  context 
of  the  linear  magnetic  induction  equation  (1.5).  The  potential-based  schemes  preserve  a 
discrete  version  of  divergence.  We  prove  numerical  stability  for  certain  versions  of  these 
schemes.  Numerical  experiments  illustrating  the  generality  of  the  approach  and  its  com¬ 
putational  efficiency  are  presented. 

The  rest  of  this  paper  is  organized  as  follows:  in  section  2,  we  describe  the  general 
form  of  constraint  preserving  potential-based  schemes.  In  section  3,  the  magnetic  induc¬ 
tion  equation  (1.5)  are  considered  and  some  stability  results  presented.  Numerical  ex¬ 
periments  are  presented  in  section  4.  This  paper  is  the  first  in  a  series  of  papers  devoted 
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to  genuinely  multi-dimensional  schemes  based  on  numerical  potentials.  Subsequent  pa¬ 
pers,  [19, 20],  will  extend  this  approach  to  non-linear  conservation  laws,  including  the 
equations  of  MHD. 


2  Potential  based  Constraint  preserving  schemes 


For  simplicity  of  the  exposition,  we  start  with  the  two  dimensional  form  of  curl  advection 
(1.3): 

(«l)t+/y  =  0, 

(«2)f-/x  =  0, 

with  the  flux  f  —  f(x,y,t,u),  where  u [  / / 1 , / / 2 )  is  the  2-vector  of  unknowns  subject  to 
divergence-free  initial  condition  divu(x,y,0)  =  0. 

We  consider  a  uniform  Cartesian  mesh  with  mesh  sizes  Ax,  Ay  in  the  x-  and  y-  direc¬ 
tions  respectively.  It  consists  of  the  discrete  cells,  =  [x{_  1  ,xi+i )  x  [y-_  1  ,y^+ 1 ),  centered 

at  the  mesh  points  (Xj,r/,)  —  f  iAx,jAi/),  (/,/)  6  Z2.  To  approximate  (2.1),  we  use  standard 
discrete  averaging  and  difference  operators 


5xaij:~- 


*i+bJ  i-bJ 


V-yal]  '■ 


IJ+ 


(2.2) 


A  word  about  our  notations:  we  note  that  the  above  discrete  operators  could  be  used 
with  indexes  l,j  which  are  placed  at  the  center  or  at  the  edge  of  the  computational  cells, 
e.g.,  I  —  i  or  I  —  i+\.  In  either  case,  we  tag  the  resulting  discrete  operators  according 
to  the  center  of  their  stencil;  thus,  for  example,  employs  gridvalues  placed  on 

the  integer-indexed  edges,  and  Wi+ 1,  whereas  5yiVj  employs  the  half-integer  indexed 
centers,  w->  1 . 

A  standard  semi-discrete  finite  volume  scheme  [18,28]  for  updating  the  cell  averages 
u y  (f)  in  (2.1)  at  time  t  can  be  expressed  as  (dropping  t  for  notational  convenience) 


—  —S  =  —  —  ( Fy  —  Vy  ^ 

Ayyfy-  Ay  IV \  Vli' 


^r(u\ )a  ■ 

dV  J  /J  Ay  '  l'J  Ay  V  w+2 

|(u2),v=Tw5.T(Fi.j-Fi.;), 


(2.3) 


where  Fx  t  .  and  FlJ.  1  are  edge  centered  numerical  fluxes,  consistent  with  the  flux  /  in 

l+  2']  l']+  2 

the  x-  and  y-  directions  respectively.  Examples  for  numerical  fluxes  include  fluxes  in  the 
viscosity  form  [26]: 


rx 

i+\,) 

Fy.  ! 
z,7+ 2 


Vxfi+ 1 +  Q Xi+  1  /x  ( U 2 ) i+ 1  j. 


(2.4) 
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where  /Z/y  :=/(xz-,yy,i, uy)  and  Qx  =  Qx(x,y,t, u)  and  =  Qy(x,y,t, u)  are  suitable  nu¬ 
merical  viscosity  coefficients.  For  example,  the  first-order  Rusanov  flux  for  (2.1)  has  the 
viscosity  form  (2.4)  with  viscosity  coefficients 


Qx,  i  .  =  max 


<4rmax 


{\xy\'\xy+i\}' 


(2.5a) 


where  A ' 1  are  the  eigenvalues  of  the  corresponding  Jacobians, 


?  =  1,2. 


(2.5b) 


The  family  of  viscous  numerical  fluxes  (2.4)  will  be  shown  to  serve  as  building  blocks  for 
the  potential-based  schemes  discussed  in  section  2.1  below. 

The  standard  finite  volume  scheme  (2.3)  may  not  preserve  any  discrete  form  of  the 
divergence  constraint  (1.4)  [11],  which  in  turn  may  lead  to  numerical  instabilities.  To 
address  this  difficulty,  we  introduce  a  family  of  genuinely  multi-dimensional  schemes 
based  on  numerical  potentials  <pi+ 1  y+i ,  defined  at  vertices  xi+i  with  the  sole  require¬ 
ment  that  these  potentials  are  consistent  with  the  differential  flux, 

cp(xryrt,u,-  ■  ■  ,u)  =  f(x,y,t,u).  (2.6) 


We  now  set  the  numerical  fluxes 


Fii+\ => 
f’+i,i=Mi+y- 


The  resulting  finite  volume  scheme  written  in  terms  of  numerical  potentials  reads 


d,  s  _  1  ,  .  _  1 

dt^U2)i’j~  Ax3*^  ~  Ax 


Ay  (2  ^!+5'/+5 

(-2^+y+\+it>i+bh0~2 


(2.7) 


The  scheme  (2.7)  is  consistent  with  (2.1)  since  the  numerical  potential  is.  There  is  a  re¬ 
markably  rich  family  of  consistent  potential-based  schemes  —  a  host  of  examples  will 
be  specified  in  the  next  subsection.  They  have  a  genuinely  multi-dimensional  structure, 
due  to  the  vertex-centered  numerical  potentials  which  include  information  in  both  nor¬ 
mal  as  well  as  transverse  directions.  Observe  that  the  potential-based  scheme  need  not 
involve  any  staggering  of  meshes.  But  before  turning  to  specific  examples  of  potential- 
based  schemes  we  describe  their  main  motivation  in  the  present  context  of  divergence- 
free  equations. 
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Lemma  2.1.  Let  uy  be  the  numerical  solution  of  the  potential-based  scheme  (2.7), (2.6).  Then, 
their  discrete  divergence,  div*,  given  by 


div*  (uy  ) 


A  +  ^yhx^y(u2)i,j/ 


(2.8a) 


is  preserved  in  time: 

^div*(uy)  =  0,  Vi,;.  (2.8b) 

Verification  of  (2.8b)  is  straightforward:  since  the  difference  operators  and  the 
averaging  operators  yx,yy  commute  with  each  other,  applying  the  discrete  divergence 
operator  div*  to  the  numerical  solution  of  (2.7)  we  find 

d  1 

^div*  (u y  )  =  -  PxSySxPy)^  =  0. 

Remark  2.1.  One  approach  in  designing  constraint  preserving  schemes  is  to  satisfy  that 
constraint  approximately :  for  example,  the  discrete  statement  of  (2.8a)  could  be  interpreted 
as  a  second-order  approximation  of  the  differential  divergence, 

div*(u/7y)  =  divu(x;,yy)  +  0(Ax2+ Ay2). 

This,  however,  requires  the  smoothness  of  the  underlying  solution.  Instead,  a  key  feature 
of  constraint  preserving  schemes  based  on  numerical  potentials  is  that  they  satisfy  exactly 
a  discrete  constraint,  so  that  their  numerical  solution  remains  on  a  discrete  sub-manifold, 
independent  of  the  underlying  smoothness. 

Lemma  2.1  shows  that  the  class  of  potential-based  schemes  satisfies  a  precise  discrete 
analogue  of  the  divergence  constraint  (1.4).  We  emphasize  that  it  applies  to  any  consistent 
numerical  potential.  Special  cases  of  the  discrete  divergence  operator  div*  (2.8a)  were 
considered  in  [29,30].  This  level  of  generality  of  the  potential-based  approach  offers  a 
major  advantage  over  earlier  studies,  and  is  explored  next. 


2.1  The  family  of  potential-based  schemes  is  rich. 

Numerical  potentials  <p  in  (2.7)  can  be  chosen  in  many  different  ways.  The  examples 
below  illustrate  the  generality  of  our  potential-based  approach. 


2.1.1  Symmetric  potential. 

A  consistent  choice  of  potential  (p  is  obtained  by  averaging  neighboring  edge  centered 
fluxes,  i.e. 


(2.9) 
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where  Fx  t  .,Fy .  A_  are  consistent  numerical  fluxes.  Higher  order  fluxes  can  be  readily 
used. 

An  explicit  computation  of  the  scheme  (2.7)  with  the  symmetric  potentials  (2.9)  leads 
to  the  revealing  form. 


d_ 

dt 

d_ 

dt 


—  4 (HxFij+i  HxFij-1 )  4At/  (^y(FxFi+1/2/j  +  FxFi_1/2/j)y 


4A  y 


(w2 )i,j  —  4Ax  (jiyFf+ij  (^(i^y^+1/2  +^^-1/2))  • 


(2.10) 


The  above  form  suggests  that  the  potential  based  scheme  (2.7)  introduces  a  special  trans¬ 
verse  correction  (by  averaging  normal  fluxes  in  the  transverse  direction)  to  the  standard 
finite  volume  scheme  (2.3).  The  above  form  brings  out  the  contrast  between  the  standard 
finite  volume  scheme  (2.3)  and  the  potential  based  scheme  (2.7)  quite  sharply. 


2.1.2  Staggered  symmetric  potential. 

A  different  consistent  potential  can  be  defined  as 

-  \ ( F* (FyV^yUi+1,J-+  \)+Fy (Vx*i+y,Vx*i+y+1 ) ),  (2-n) 

where  Fx,Fy  are  consistent  two-point  numerical  fluxes.  This  approach  is  equivalent  to 
the  symmetric  potential  (2.9)  for  linear  equations  with  constant  coefficients.  However,  it 
leads  to  a  different  scheme  for  equations  with  variable  coefficients  and  for  equations  with 
non-linear  fluxes. 


2.1.3  Diagonal  Potential. 

A  completely  different  form  of  the  potential  is  defined  by 

Qi+\,j+\  =  \  (f*K/Vum,;+ 1 )  +  Fy(uy,u/+1,;-+1 )),  (2.12) 

where  Fx  and  Fy  are  consistent  two-point  numerical  fluxes.  It  is  straightforward  to  extend 
this  form  for  any  2fc-point  numerical  fluxes.  This  form  of  the  potential  is  isotropic  and 
leads  to  a  compact  form  of  the  scheme  (2.7). 

2.1.4  Mixed  Potential. 

A  slightly  different  form  of  the  diagonal  potential  is  obtained  as, 

<t>i+y+\= 

—  ^  (-F  (uJ,/'/UJ+l,/'+l  )  +  F  (Ui,j+1  rUi+l,j  )  +  f  ^ ( U ! + 1 ,/' / U !,/' + 1  )  +  f ( U !,/' / U i+ 1 ,/' + 1  X^,13) 

where  Fx,Fy  are  consistent  two-point  numerical  fluxes. 
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2.2  Second  order  potential-based  schemes 

The  order  of  accuracy  of  the  scheme  (2.7)  is  related  to  the  choice  of  numerical  fluxes 
Fx ,FV  used  to  define  the  potentials.  The  discrete  divergence  operator  div*  (2.8a)  and  the 
difference  and  averaging  operators  (2.2)  are  second-order  accurate.  Therefore,  overall 
second-order  spatial  accuracy  is  obtained  by  using  standard  non-oscillatory  piecewise 
linear  reconstructions  in  each  cell.  We  follow  the  reconstruction  procedure  proposed 
in  [17]. 


2.2.1  Second  order  non-oscillatory  reconstruction. 


The  cell  averages  uy  are  used  to  define  the  piecewise  bilinear  reconstruction. 


(2.14a) 


The  numerical  derivatives  in  x—  and  y—  directions,  denoted  respectively  by  prime  and 
backprime  are  given  by  the  standard  limiter. 


U  y  =  minmod(uj+ij  -  Uy,0.5(um,;-  -  -  u /-y), 

u \j  =  minmod(uy+i  -  Uy,0.5(Uy+i  -  Uy_i ),uy  -  Uy_i ) . 
The  minmod  function,  defined  as 


(2.14b) 


minmod  ( a, b,c) 


sgn(a)min{\a\,\b\,\c\}, 

0, 


if  sgn(a)  —sgn(b)  —  sgn(c), 

otherwise. 


is  a  standard  van-Leer  limiter;  other  standard  limiters  can  also  be  used  in  this  context. 
The  use  of  limiters  ensures  that  the  reconstruction  of  each  unknown  is  non-oscillatory 
e.g.,  total-variation  diminishing,  consult  [28]  and  the  references  therein.  In  the  sequel,  we 
will  also  need  the  following  reconstructed  corner  pointvalues, 

utj  =  V  i,j{xi+hyj),  uJJ  =  pij  ,y}), 

. 2  2  (2.14c) 

=  Vi,i  (xuyj+\ )’  ui,j  =  Vi.)  (xuVj- 1 )  • 

Given  the  any  consistent  two-point  fluxes  F,G  ,  a  second  order  flux  based  on  a  midpoint 
rule  to  compute  edge  integrals  takes  the  from. 


FhrFX«>'u 


F  u/+i,y ) 


2 


7+i '  * 


(2.14d) 


The  above  fluxes  are  used  to  define  the  potentials  in  the  scheme  (2.7)  resulting  in  an 
overall  second  order  accurate  discretization  of  (2.1). 


Remark  2.2.  In  order  to  achieve  third  and  even  higher  order  accuracy  in  space,  we  need  to 
redefine  the  averaging  and  difference  operators  given  in  (2.7)  to  higher  than  second  order 
accuracy.  The  potentials  can  then  be  expressed  in  form  of  fluxes,  based  on  third  and  even 
higher  order  (W)ENO  type  reconstructions,  [15,25].  This  program  for  designing  schemes 
of  arbitrary  orders  of  accuracy  will  be  considered  in  a  forthcoming  paper. 
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2.3  A  divergence  preserving  viscous  discretization. 

The  two  dimensional  form  of  the  divergence  preserving  viscous  equation  (1.7)  is 

(  (wi)f  H“  fy  —  {u{)yy  —  (u,2)xy/ 

1(^2 )t~  fx~  {ul)xx  —  {U\ )xy- 


(2.15) 


Note  that  the  viscosity  in  (2.15)  is  of  the  curl  (curl)  type. 

We  combine  the  potential  based  discretization  of  the  flux  terms  in  (2.1)  with  a  simple 
genuinely  multi-dimensional  discretization  of  the  viscous  term  to  obtain  a  divergence 
preserving  scheme  for  (2.15).  Below,  we  employ  the  standard  notations  for  forward, 
backward  and  centered  divided  differences. 


j^±  _^z(ai±l,j  ai,j) 

x  a,'j  ~  Ax 


n±  _  (ai,j±~L  ai,j ) 

uvai'i-  a y 


D°  n  -  — _ — 

Lyx,yui,]  ~ 


^x,yai,j  +  ^x,yai,i 


The  divergence  preserving  scheme  for  (2.1)  is 

dl  1 

)i,j  —  —  -^yV'xtyi,)  —  DxDy(u2)i,j 

+  7  (DyDy  (Wl)^'+l ,/+2Dy  Dy  (wl )i,j  +  DyDy  (t/i)/-i,yY 

d  i  (2-16) 

^(^2 )i,j  ~  —  DxDy  (ul)i,j 

+ 1  (u2)i,j+i  +2D+ D~  (u2)i,j  +  D+ D~  (u2)i,j-i  )  • 

The  above  scheme  is  a  consistent  discretization  of  (2.15).  A  straightforward  calcula¬ 
tion,  together  with  Lemma  2.1  shows  that  the  scheme  (2.16)  preserves  the  discrete  di¬ 
vergence  operator  div*.  Note  that  the  viscous  terms  are  discretized  in  a  genuinely  multi¬ 
dimensional  manner  in  (2.16). 


2.4  Curl  preserving  discretization 

The  g rad  advection  equation  (1.8)  in  two  space  dimensions  with  flux  f  =  f(x,y,t,\x)  is 
given  by 

J(ki),+A=0,  (2.17) 

\(u2)t+fy  —  0. 

We  follow  the  strategy  of  the  previous  sections,  seeking  a  consistent  vertex  centered 
potential,  <pi+ 1  y+i  •  The  corresponding  potential-based  scheme  reads 


(2.18) 
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The  scheme  preserves  the  following  discrete  curl  operator, 
d  11 

^curl*  (u  i,j)  =  0,  curl*  (uy )  =  —  /^(wOy  -  — ^(w2)y  •  (2.19) 

The  proof  follows  along  the  lines  of  lemma  (2.1).  The  potentials  can  be  defined  in  a  man¬ 
ner,  analogous  to  the  divergence  preserving  scheme  (2.7).  Indeed,  the  two-dimensional 
divergence  and  curl  preserving  equations,  (2.1)  and  (2.17),  and  their  corresponding  potential- 
based  schemes,  (2.7)  and  (2.18),  are  dual  to  each  other. 

2.5  Divergence  preserving  schemes  in  three  dimensions 

The  three  dimensional  divergence  preserving  equations  (1.3)  with  flux  vector  f—  (/i,/2,/3), 
are  explicitly  written  as, 

(ul)t  +  (/3 )y  ~  (fl)z  =  0, 

(u2)t  +  (fi)z-(f3)x  =  0,  (2-20) 

(u3 )t  +  (fl)x  ~  (/l)y  =  0, 

It  is  straightforward  to  extend  the  potential  based  framework  of  this  section  and  design  a 
divergence  preserving  scheme  for  the  above  equation.  Let  Sx,Sy  and  5Z,  px,py  and  yz  de¬ 
note  the  difference  and  average  operators  in  the  x,y  and  z  directions  respectively.  Define 
a  uniform  grid  in  all  three  directions  (x/,yy,z^)  =  (zAx,;Ay,/cAz)  with  mesh  sizes  Ax,A y 
and  Az.  Also,  denote  the  cell  C^jc  —  [xf_i,x/+i)  x  [y-_i,y.+  i)  x  [zk_i,zk+i)  and  the  cell 
average  of  the  unknown  over  the  cell  as  u^ . 

We  need  to  define  three  vertex  centered  potentials  (</>/) -+ 1  -+  ik+i  with  l  — 1,2,3,  such 
that  they  are  consistent ,  i.e, 

(<pi)i+y+ty+i(x,y,zj,ur--,u)=fi(u),  l  =  1,2,3.  (2.21) 

The  divergence  preserving  scheme  in  three  dimensions  is  defined  in  terms  of  the  poten¬ 
tials  as 

(j!  i 

d  \  1 

)i,j,k  ~  ~~^.^z'H-x'H-y((pl)i/j/k  +  -^>xV'yV'z{<p3)i,jkt  (2.22) 

d  i  | 

3)i,)X  —  ~  ~^^xyyyz{^2)i,j/kJr-^^y}ixyz{(pl)i/j/k' 

Arguing  along  the  lines  of  lemma  2.1  we  now  state  the  following  divergence  preserv¬ 
ing  property. 

Lemma  2.2.  Let  u ^  be  the  numerical  solution  of  the  potential-based  scheme  (2.22),(2.21).  Then , 
their  discrete  divergence ,  div*,  gzizen  by 

1  1  1 

div  (u^ )  =  y y y W 1 ) /ry,/c  T-  Pxhz^yf^'l) i,j,k  PyPx$z{u 3 ) /,/,/c / 


(2.23a) 
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is  preserved  in  time: 

jjdiv*  (u«# )  =  Vi,;,  A:.  (2.23b) 

There  is  a  rich(-er)  class  of  3D  consistent  potentials  which  can  be  defined  in  a  manner 
analogous  to  the  two  dimensional  case.  As  an  example,  we  define  the  three  dimensional 
form  of  the  symmetric  potential  (2.9),  (p\  as  follows. 


g  ((^1  )i,j+\jc  +  i^l)i+l,j+\,k  +  (^1  )«',;+ 5^+1  “*■  (^1  )«+!,;+ 5^+1 ) 


+  ■ 


*?Wi  +  (f?) 


-i  +  tf 


z+l,/,A:+ 1  '  vrl  ljc+ \  '  Vrl  /i+l,/+l/fc+ 1 


i  +  W 


where  F\,F\  are  numerical  fluxes  consistent  with  the  flux  f\  in  the  y-  and  z-  directions 
respectively.  The  potentials  (pi,(p3  can  be  similarly  defined. 


2.6  Time  stepping 

The  constraint  preserving  schemes  discussed  so  far  are  semi-discrete  schemes  and  need 
to  be  coupled  with  suitable  time  integration  routines.  Consider,  for  example,  the  2D 
potential-based  scheme  (2.7)  at  time  level  t  —  tn, 


— u  -(tn)  —  E”.  E — 


1 


(2.24) 


The  simplest  time  stepping  is  the  first-order  accurate  forward  Euler  scheme. 


un+x  =u?-  —  E^- 
uz,/  uz,7  Ez,/' 


(2.25) 


where  the  time  step  A  tn  is  limited  by  a  suitable  CFL  condition.  Second-order  temporal 
accuracy  can  be  obtained  using  a  SSP  Runge-Kutta  method  [14], 


—,,n  A  +nVn 

U/,y  -  Ui,j  -  M 


>)<>)<  5)« 

ui,j  =+/- 


-  Af"E"y, 


u 


n+1 _ 

V  ' 


(2.26) 


An  alternative  first-order  accurate  genuinely  multi-dimensional  time  stepping  is  the  ex¬ 
tended  Lax-Friedrichs  type  time  stepping. 


u 


n+1 _ 

y  i 


:  (4u£/  +  um,j  +  ulj+ 1  +  uf-i,;  +  ulj- 1 


i  —  Af"E"y. 


(2.27) 


13 


3  Schemes  for  the  Magnetic  induction  equation 


The  preceding  description  on  constraint  preserving  schemes  is  very  general.  In  order  to 
provide  some  concrete  stability  estimates  and  perform  numerical  experiments,  we  focus 
on  the  two-dimensional  form  of  the  magnetic  induction  equations: 

f  («l)f  +  (V2U1  -  VlU2)y  =  0, 

\(U2)t-(V2Ul-ViU2)x  =  0, 


with  the  magnetic  field  u  =  (ui,u2)  and  a  given  velocity  field  v  —  {v\,v2). 

In  order  to  complete  the  divergence  preserving  potential  based  scheme  (2.7),  we  need 

to  specify  numerical  fluxes  Fx  L  .,Fy .  L.  The  simplest  available  two  point  flux  is  the  aver- 

l+2'J  V+2 

age  flux: 


rx 


f bj  f i+i,y  -pV  _ fhj+Aj+i 

7  tj+\ 


(3.2) 


2  7  ii+h  2 

where  /( u)  =  V2U\  —  V\U2-  Using  the  above  flux,  together  with  the  symmetric  potential 
(2.9)  results  in  the  potential-based  scheme: 


1 


dt 

d 

Tt 


(«2k 


1 

';_2Ax 


(3.3a) 


where 


fi,j  Hxfi,j  =  |  {fi+bj+Zfij+fi-lj),  fi,j  Fyfi,j  =  |  (/y+1  +2/y  +fi,j-l)  •  (3.3b) 

The  scheme  (3.3a),  which  coincides  with  the  symmetric  scheme  proposed  in  [30],  preserves 
the  discrete  divergence, 

div*  (uy(f))  =div*  (uy  (0)),  div*(uy)  :=-F}iy5x(u1)i/j  +  ^F}ix5y(u2)i,j.  (3.4) 


We  will  show  that  the  potential-based  symmetric  scheme  (3.3)  is  L2-stable.  Our  proof 
highlights  the  role  of  (discrete)  divergence-preserving.  To  motivate  the  numerical  sta¬ 
bility,  we  first  provide  the  corresponding  L2  well-posedness  statement  in  the  continuous 
case. 


Lemma  3.1.  [1 1  ].  Let  u  be  the  weak  solution  of  the  magnetic  induction  equations  (3.1)  subject  to 
divergence  free  initial  data,  uo6  L2(R2),  and  a  convective  velocity  field  v  (vi,v2)  E  C1(R+,R2). 
Then,  u  satisfies  the  apriori  energy  boundf 


d_ 

dt 


vllci 


(3.5) 


+ We  use  X Y  to  denote  the  estimate  X  <  CY,  where  C  is  a  constant  which  may  depend  on  t  but  otherwise  is 
independent  of  the  solution.  Ax, Ay  etc. 
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To  verify  (3.5)  we  follow  [11].  Adding  to  (3.1)  a  zero  source  term  which  is  proportional 
to  vanishing  divergence,  we  obtain. 


(«l)f  +  (v2Ul  -  VlU2)y  =  -Vl(Ul)x  -  Vl(u2)y, 
(U2)t  ~  {V2U\  ~  ViU2)x  =  ~V2(U1)X  -  V2(u2)y, 


(3-6) 


so  that  after  straightforward  simplifications,  (3.1)  recast  into  the  symmetric  form, 

{Ul)t  +  Vl{Ul)x  +  V2(Ul)y=-(V2)yUl  +  (Vl)yU2, 

(u2)t  +  Vl(u2)x  +  V2(u2)y  =  (v2)xUl  -  (vi)xu2. 

The  desired  L2-bound  follows  by  applying  the  energy  method  for  the  symmetric  form  of 
the  equations  in  (3.7)  with  bounded  low-order  term,  ||v||ci  <  o°. 

We  now  turn  to  show  that  the  potential-based  symmetric  scheme  (3.3a)  satisfies  a 
discrete  version  of  the  L2  energy  estimate  (3.5). 

Lemma  3.2.  Let  u ^  be  the  solution  of  the  semi-discrete  potential  based  symmetric  scheme  (3.3) 
subject  to  divergence-free  initial  data  div*(u(jq,i/y,0)=0  and  velocity  field  veC1(R+,R2).  Then, 
the  following  L1 -energy  estimate  holds, 

4|MIl2~  IMIc^HIl2'  ||u||22  :=A*Ay£(Mi)?.  +  (tt2)?..  (3.8) 

t  A  A  A  .  . 


Proof.  We  mimic  the  proof  of  the  continuous  estimate  (3.5).  We  begin  by  writing  the 
discrete  divergence  div*  in  the  form 


where  ux  and  0  are  the  averages  (3.3b).  The  discrete  divergence  preservation  (3.4)  tells 
us  that  div*  (u hf)  =  0:  adding  a  multiple  of  this  vanishing  divergence  as  a  discrete  source 
term  to  (3.3a)  yields  a  discrete  analogue  of  (3.6)), 


-D %  -  MjD^uDij  -  {vfijDyUfi;, 


dt 

d 


(«2)y  D°xff. -  (p2)i,jDj(«i)i,j  -  0 *2)y DS(fif)y • 


^0 1  ,-.x  \ 


(3.9) 


dt 


Substituting  the  explicit  form  oi  f—v2U\  —  u2V\  in  (3.9),  one  obtains  after  straightforward 
manipulations,  the  following  discrete  version  of  the  symmetric  form  of  the  equations  in 
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(3.7), 


./  =  -(»i)«DS(“i)«  _  ('’DijKlij+i +D,  (*2)ij(“i)i,/-i 

+  j(Dy+  (Sf  )y  (fil)«+l  +  Dy“  W)y  (S|)y-i), 

(w2 )i,j  ~  ~  ~  (v2 )i,jDy(U2)i/j  +  -  (wfyhj  (™l)i+l,j  + 

~  2  (^1  )l/ (^2)*,/+!  (v\)i,j(U2)i,j- 


d 

dt 


(3.10) 


We  conclude  with  energy  method:  summing  by  parts  the  first  equation  in  (3.10)  against 
AxAy(ui)ifj  and  the  second  equation  against  AxA  1/(^2)/,/  and  using  Cauchy's  inequality, 
we  obtain  the  L2  energy  estimate  (3.8).  □ 


Note  that  preserving  the  discrete  divergence  (3.4)  plays  a  crucial  role  in  the  Instability 
of  the  potential-based  scheme  (3.3).  Since  the  scheme  is  based  on  centered  stencil  it  is  un¬ 
conditionally  unstable,  when  combined  with  the  forward  Euler  or  second-order  Runge 
Kutta  (RK)  time  stepping  (2.25),  (2.26).  Stability  can  be  achieved  by  using  higher-order 
(>  3)  RK  time  stepping,  e.g.,  [27].  Alternatively,  a  standard  way  to  stabilize  centered- 
based  schemes  is  achieved  by  adding  numerical  diffusion  (2.4).  A  simple  Rusanov  type 
numerical  diffusion  operator,  (2.5),  is  used  to  modify  the  numerical  fluxes  (3.2),  yielding 
the  viscous  numerical  fluxes, 

1 

Fi+y  =  =  -  (fi,j+fi+ y)  +  max{  |  (i>i)y  |,  |  (oi)i+i,;  | }  {(u2)i+1/j  -  (w2)y) 

\ 

F(j+i  =Fy(ui,;,ui+i/;)  =  ^(hj+hj+i) -max{\(v2)ifj \r\(v2)i,j+i 1}  ((wi)y+i  -  («i)y). 

(3.11) 

The  diffusive  terms  involve  the  maximum  wave  speeds  in  each  direction.  Other  diffu¬ 
sion  operators,  like  the  standard  upwind  diffusion,  can  also  be  used.  Once  the  diffusive 
numerical  fluxes  are  set,  one  can  define  the  corresponding  numerical  potential  and  com¬ 
plete  the  potential-based  scheme  (2.7).  We  are  unable,  however,  to  prove  that  the  scheme 
(2.7)  with  numerical  diffusion  (like  in  (3.11))  is  L2-energy  stable.  Nevertheless,  numer¬ 
ical  experiments  in  the  sequel  suggest  that  the  potential-based  scheme  with  numerical 
diffusion  is  energy  stable. 


4  Numerical  Experiments 

We  test  the  constraint  preserving  schemes  for  the  magnetic  induction  equation  (3.1)  in 
this  section.  The  following  four  schemes  are  considered: 

All  schemes  are  updated  in  time  with  a  CFL  number  of  0.45. 
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RUS  Standard  first-order  Rusanov  scheme,  (2.3),(3.11). 

CPR  Constraint  preserving  scheme  (2.7)  with  Rusanov  flux. 

CPR2  Second-order  (both  space  and  time)  constraint  preserving  scheme,  (2.14),  based  on  Rusanov  flux  (2.5). 
CPS  Constraint  preserving  symmetric  scheme  (3.3a)  with  the  Runge-Kutta  time  stepping. 

4.1  Numerical  experiment  #1:  rotating  hump. 

This  test  case  is  a  benchmark  for  testing  schemes  for  multi-dimensional  advection  [11,30]. 

We  consider  the  two  dimensional  magnetic  induction  equation  (3.1)  with  the  velocity 
field  (zq,r>2)  =  (— y,x).  The  exact  solution  can  be  calculated  as 

u(x,y,t)  =  R(t)u0(R(-t)(x,y)),  (4.1) 

where  R(t)  is  a  rotation  matrix  with  angle  t. 

We  consider  the  divergence  free  initial  data: 

U0 (X,y)  —4  (x ~yi)  e-H(x-\)2+y2) ,  (4.2) 

and  the  computational  domain  [—1,1]  x  [—1,1].  The  exact  solution  (4.1)  is  a  smooth  hump 
(centered  at  (^,0))  rotating  about  the  origin  and  completing  one  rotation  in  time  t  —  2n. 

Non-reflecting  Neumann  type  boundary  conditions  are  used. 

The  approximate  solutions  at  time  t  =  2n  on  a  100  x  100  mesh  are  shown  in  figure  1. 

We  show  the  norm  ||  u  ||  =  yju\  +  u\  with  four  different  schemes.  The  constraint  preserving 
schemes  are  based  on  the  symmetric  potential  (2.9).  The  figure  shows  that  the  standard 
RUS  scheme  does  a  poor  job  of  approximating  the  rotating  hump.  The  magnitude  of 
the  hump  is  smeared  considerably  and  the  shape  is  distorted.  Unphysical  waves  are  also 
generated.  In  sharp  contrast,  the  constraint  preserving  schemes  approximate  the  solution 
quite  well.  The  first-order  CPR  scheme  smears  the  solution  somewhat  (note  different 
scales  in  the  figures)  but  the  shape  is  still  maintained.  The  second-order  CPR2  and  the 
CPS  schemes  resolve  the  solution  very  sharply.  The  smearing  is  reduced  considerably 
and  the  shape  is  maintained.  The  results  suggest  a  strong  connection  between  divergence 
preservation  and  the  numerical  performance.  This  link  is  quantified  in  Table  1  where  we 
tabulate  the  discrete  divergence  div*  (2.8a)  in  L2,  generated  with  all  the  four  schemes  at 
time  t  —  2 re.  Note  that  since  the  initial  data  is  divergence  free,  the  exact  divergence  will 
be  zero  at  all  times. 

Table  1  shows  that  the  standard  RUS  scheme  generates  divergence  errors  of  the  or¬ 
der  of  the  truncation  error  on  coarse  meshes.  However,  the  scheme  is  unstable  on  fine 
meshes  and  crashes  on  a  400  x  400  mesh.  The  blow  up  of  RUS  scheme  based  on  400  x  400 
mesh  points  was  preceded  by  a  large  increase  in  the  divergence,  indicating  a  possible  con¬ 
nection  between  the  two.  The  constraint  preserving  schemes  result  in  very  small  diver¬ 
gence  errors,  mostly  due  to  boundary  effects  (no  special  divergence  cleaning  is  applied 
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(c)  CPR2  (d)  CPS 

Figure  1:  ||u||  at  time  t—2n  in  the  rotating  hump  experiment  with  four  different  schemes  on  a  100x100  mesh. 


Grid  size 

RUS 

CPR 

CPR2 

CPS 

50/50 

2.15e-l 

6.0  7e-7 

1.27  e— 7 

7.0e  —  7 

100/100 

8.15e-2 

3.04e-9 

2.9e  —  8 

1.5e-10 

200/200 

2.2e+3 

5.34e-12 

2.88e  — 14 

2.8e-U 

400/400 

diverges 

5.53e-14 

7.1e  — 15 

7.8e  — 15 

Table  1:  Absolute  errors  in  L2  for  ||div*||  at  time  t  —  2n  for  the  rotating  hump  with  four  different  schemes. 
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at  the  boundary).  The  errors  on  coarse  meshes  are  very  small  and  converge  to  zero  quite 
rapidly  As  expected,  the  errors  are  of  the  order  of  machine  precision  on  fine  meshes. 

The  operator  div*  is  a  particular  discrete  form  of  the  divergence  operator.  We  have 
shown  that  the  constraint  preserving  scheme  (2.7)  preserves  this  particular  form  of  diver¬ 
gence.  A  natural  question  is  what  happens  when  a  different  form  of  discrete  divergence 
is  used.  We  consider  the  standard  centered  discrete  divergence  operator. 


div(u itj)  -  D°(wi)y  +  D°(u2)i,j. 


A  simple  calculation  shows  that  div  and  div*  differ  by  0(Ax2+ Ay2).  Therefore,  pre¬ 
serving  div*  would  only  imply  that  errors  in  div  behave  like  the  square  of  the  truncation 
error.  This  issue  is  explored  quantitatively  and  the  results  are  shown  in  Table  2.  The 


Grid  size 

RUS 

CPR 

CPR2 

CPS 

50/50 

2.16e-l 

1.48e-2 

1.43c  — 2 

1.8c  — 2 

100/100 

8.29e-2 

3.2e-3 

1.32c -3 

1.6c  — 3 

200/200 

1.99e+3 

6.26e-4 

6.7c— 5 

2.0c -4 

400/400 

blow-up 

1.67e-4 

4.4c -5 

5.6c -6 

Table  2:  Absolute  errors  in  L2  for  1 1 div 1 1  at  time  t  =  2n  for  the  rotating  hump  with  four  different  schemes. 

table  shows  that  the  errors  in  the  standard  divergence  div,  generated  by  the  constraint 
preserving  schemes,  are  low  and  are  consistently  lower  than  the  expected  square  of  the 
truncation  error.  Consequently,  we  conclude  that  using  a  divergence  preserving  scheme 
will  lead  to  lower  errors  for  other  discrete  forms  of  divergence. 

In  the  above  discussion,  the  divergence  preserving  schemes  were  based  on  the  sym¬ 
metric  potential  (2.9).  We  use  all  the  four  potentials  described  in  section  2  with  the 
first-order  CPR  scheme  and  show  ||u||  for  the  approximate  solution  at  time  t  —  tt/4  on 
a  100  x  100  mesh  in  figure  2.  The  figure  clearly  shows  that  different  choices  of  poten¬ 
tial  resulted  in  very  similar  numerical  approximations.  This  was  also  seen  in  other  ex¬ 
periments,  indicating  the  robustness  of  our  approach  with  respect  to  varying  choices  of 
potentials.  However,  there  were  some  boundary  instabilities  for  the  diagonal  potential 
when  long  time  scales  were  considered.  This  fact  requires  careful  investigation  in  the 
future. 


4.2  Numerical  experiment  #2:  discontinuous  test  case. 


The  rotating  hump  involved  smooth  solutions.  However,  we  can  expect  discontinuous 
solutions  (particularly  in  MHD  models).  We  test  the  constraint  preserving  scheme  on  a 
numerical  experiment  involving  discontinuities  in  the  solution.  The  initial  data  is  given 
by,  [11], 


ui(x,y)=i4(x,y)  = 


2 

0 


if  x  >  y, 
otherwise. 
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(c)  Diagonal  (d)  Mixed 

Figure  2:  Rotating  Hump:  ||u||  at  time  t  =  zr/4  on  a  100x100  mesh  with  the  CPR  scheme  and  different 
potentials. 
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and  the  velocity  field  is  a  constant  v=  (1,2).  The  exact  solution  is  the  initial  discontinuity 
moving  along  the  diagonal, 

VL(x,y/t)  =  u°(x  —  t/y—2t). 

We  consider  the  computational  domain  [—2,2]  x  [—2,2]  and  use  extrapolated  Neumann 
boundary  conditions.  A  one  dimensional  slice  (at  y  =  0.0)  of  the  solution  component  U\ 
computed  with  the  RUS,  CPR  and  CPR2  schemes  at  time  t  —  0.5  on  a  100  x  100  mesh  is 
shown  in  figure  3.  The  figure  shows  that  the  standard  RUS  scheme  leads  to  a  large  over- 


Figure  3:  Numerical  experiment  2:  ui(x,0,0.5)  on  a  100x100  mesh  with  different  schemes. 

shoot.  This  fact  was  first  observed  in  [11].  Furthermore,  the  discontinuity  is  smeared 
to  a  large  extent.  The  CPR  scheme  reduces  the  overshoot  considerably  and  resolves 
the  discontinuity  with  very  little  smearing.  However  there  are  small  amplitude  oscil¬ 
lations,  showing  that  the  constraint  preserving  scheme  is  not  necessarily  total-variation 
diminishing  ( TVD ),  although  the  exact  solution  in  this  particular  experiment  is  TVD. 
The  second-order  CPR2  scheme  resolves  the  discontinuity  more  sharply  and  the  oscilla¬ 
tions  are  also  reduced.  The  results  show  that  the  constraint  preserving  schemes  are  not 
diffusive  enough  in  this  case.  A  simple  method  for  increasing  the  diffusion  without  af¬ 
fecting  the  constraint  preserving  abilities  is  to  use  the  genuinely  multi-dimensional  Lax- 
Friedrichs  type  time-stepping  (2.27)  with  the  CPR  scheme.  We  term  this  scheme  and  its 
second-order  (in  space)  version  as  aCPR  and  aCPR2  scheme  respectively.  The  results  are 
shown  in  figure  4.  The  aCPR  scheme  removes  the  oscillations,  at  the  cost  of  smearing  the 
discontinuity.  The  spatially  second-order  aCPR2  scheme  captures  the  discontinuity  more 
sharply  and  without  any  noticeable  oscillations.  This  alternative  time  stepping  provides 
a  simple  recipe  of  modifying  the  constraint  preserving  schemes  to  eliminate  unphysical 
oscillations. 
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Figure  4:  Numerical  experiment  2:  ui(x,0,0.5)  on  a  100x100  mesh  with  different  schemes. 

Remark  4.1.  The  potential  based  scheme  (2.7)  is  slightly  more  expensive  than  its  build¬ 
ing  block,  the  standard  finite  volume  scheme  (2.3).  However,  the  overall  cost  is  still  quite 
low.  The  simplicity  and  generality  of  this  approach  renders  it  considerably  cheaper  and 
easier  to  implement  than  competing  constraint  preserving  frameworks.  The  extra  com¬ 
putational  cost  is  justified  by  the  considerable  increase  in  stability  and  resolution. 


5  Conclusion 

Evolution  equations  (1.1)  with  an  intrinsic  constraint  (1.2)  are  considered.  Examples  in¬ 
clude  the  divergence  preserving  equations  (1.3)  and  the  curl  preserving  equations  (1.8). 
Standard  finite  volume  schemes  (2.3)  do  not  necessarily  preserve  discrete  versions  of  the 
constraint  and  may  be  unstable. 

We  design  finite  volume  schemes  for  (1.1)  that  preserve  discrete  forms  of  the  con¬ 
straint  (1.2).  The  schemes  are  based  on  vertex  centered  numerical  potentials.  The  resulting 
scheme  is  genuinely  multi-dimensional  and  constraint  preserving.  The  class  of  potential 
based  schemes  is  very  rich.  Potential  based  schemes  are  presented  for  both  the  diver¬ 
gence  preserving  equation  (1.3)  and  curl  preserving  equation  (1.8).  Constraint  preserving 
schemes  for  the  equations  with  viscosity  (1.7)  are  also  proposed.  Second-order  accuracy 
is  obtained  by  employing  non-oscillatory  piecewise  polynomial  reconstructions. 

The  magnetic  induction  equations  in  two  dimensions  (3.1)  are  considered  in  detail. 
A  divergence  preserving  potential  based  scheme  is  shown  to  be  L2-stable.  Numerical 
experiments  demonstrating  the  robustness  and  computational  efficiency  of  the  constraint 
preserving  schemes  are  presented.  They  show  that  the  schemes  perform  considerably 
better  than  standard  schemes. 
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The  main  advantage  of  this  new  approach  is  its  simplicity  and  generality.  This  pa¬ 
per  is  the  first  in  a  series.  Subsequent  papers  include  ones  describing  potential  based 
genuinely  multi-dimensional  schemes  for  systems  of  conservation  laws  [19]  and  diver¬ 
gence  preserving  schemes  for  the  ideal  MHD  equations  [20].  Extending  the  potential 
based  schemes  to  higher  than  second  order  accuracy  and  to  unstructured  grids  will  be 
considered  in  future  papers. 
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